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We study the chemical potential of water as a function of charge based on perturbation theory. 
By calculating the electrostatic-energy fluctuations of two states (fully charged and uncharged) we 
are able to determine accurate values for the dependence of the chemical potential on charge. We 
find identical results for the chemical-potential difference of fully charged and uncharged water 
from overlapping-histogram and acceptance-ratio methods and by smoothly connecting the curves 
of direct exponential averages. Our results agree with those of Rick and Berne (J. Am. Chem. 
Soc., 1994 , 116, 3949) with respect to both the chemical-potential difference and its dependence on 
the charge coupling parameter. We observe significant deviations from simple Gaussian-fluctuation 
statistics. The dependence on the coupling parameter is not quadratic, as would be inferred from 
linear continuum models of electrostatics. 


I. INTRODUCTION 

The accurate calculation of thermal properties of fluid systems such as the free energy or the entropy using computer 
simulations poses serious difficulties. Unlike mechanical quantities (e.g., energy, pressure) which can be expressed as 
averages of functions of the phase-space .coordinates, thermal quantities are related to the partition function and thus 
to the volume of accessible phase space.El 

We recently pursued an approach based on perturbation theory for the chemical potential which relates the free- 
energy difference of two states to fluctuations of the two equilibrium systems™ This has the advantage that along 
with free-energy data other useful information can be extracted from equilibrium computer simulations. Motivated by 
the Born-model prediction of quadratic dependence of the free energy on ionic chargep we studied the free energy of 
simple ions in water. Within the framework of perturbation theory, the Born model suggests validity of a second-order 
treatment. with the ionic charge as a coupling parameter and, correspondingly, Gaussian statistics of the electrostatic 
potential. Em Our work indeed supported this view; but it also showed that to get accurate values for the free energy 
of ionic solvation it is important to combine fluctuation information of different charge states. 

In this work we study the chemical potential of water. This quantity plays a central role in many processes of 
physical chemistry and biophysics. For instance, the chemical potential of water is the driving force in osmotic 
equilibria. To study inhomogeneous aqueous systems or mixtures such as macromolecular solutions or crystals using 
grand-canonical-ensemble methods, accurate values for the chemical potential of the particular water model are 
required. 

We will be concerned with the electrostatic contributions to the chemical potential. Efficient methods such as 
test-particle insertionD are available for the contributions related to the particle volume. We will study the process of 
uncharging a water molecule with the van der Waals interactions. unmodified. Fluctuation statistics and, for reference, 
Bennett’s overlapping-histogram and acceptance-ratio methodsQ will be used to calculate the chemical potential. We 
also use a geometrical method which smoothly connects the direct exponential averages. We compare our results for 
the chemical potential with those of Rick and BerneQ obtained from thermodynamic integration. 

We will first develop, the theoretical framework. The correction for finite system size will be discussed. As in 
our previous studies, ullj we will use Ewald lattice summation and a generalized reaction-held (GRF) method for the 
electrostatic interactions and consistently add the self-interactions as finite-size corrections. We will then describe the 
computer simulations and discuss the results. 
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II. CALCULATION OF CHEMICAL POTENTIALS 


A. Fluctuation statistics 


Various methods are available to calculate free-energy differences between different states of fluid systems using computer 
simulations (see refs |I|,jll[|l2| for reviews). Here we are concerned with determining the electrostatic contributions to the free 
energy of water. The chemical-potential difference between a charged and an uncharged water molecule in bulk water is 

calculated. Our goal is to relate A[i ex to fluctuations in the electrostatic energy of single water molecules. ri 
The potential-distribution theorem for the excess chemical potential fj, ex forms a convenient starting point B 

fj, ex (Xi) — fi ex (\o) = —fefiTln (exp {—/3[w(Ai) — m(Ao)]}} Ao . (1) 

Ao and Ai are coupling parameters representing the two charge states. (3 is the inverse temperature 1/fcsT. The thermal 
configuration-space average in the charge-state A is denoted by (.. .)a- u( A) is the charge- and configuration-dependent interac¬ 
tion energy of a water molecule carrying charges (1 — A )qo and (1 — A )qH on oxygen and hydrogen atoms, respectively. A = 0 
and 1 correspond to the fully charged and uncharged state, respectively. We will only discuss results for three-point models of 
water. Extensions to water models carrying four and more partial charges are trivial. 

The energy w(A) contains the electrostatic interactions u e i(X) with all other water molecules and a self-interaction u 3 ( A), 

m(A) = Uei{X) + m s (A) , (2) 


where 


u e i{X ) = (1 — A) [qo4>o + qn ( 4>h 1 + ■ (3) 

4>o, and 4>h 2 are the electrostatic potentials at the sites of the oxygen and hydrogen atoms. The self-interaction u s {X) 

depends on the effective electrostatic interaction <p(r) used in the computer simulations, 

Us{X) = Xfq a qp^{r al3 ) , (4) 

a,(3 


where the double sum extends over pairs of charge sites on a water molecule. VP is averaged over the solid angle fl, 


^(rais) 


/ dfl 

<p(r a/3 ) - 7^7 

Jci 

L l r IJ 


(5) 


The perturbation expression eq |l] can be used directly if the Ai and Ao states are close. Otherwise, the average, in eq [l] is 
dominated by the poorly sampled tails of the distribution of At = k(Ai) — u(Ao). As in our previous workJ3tj we use a 
perturbative expansion. Application of the cumulant expansionli3 with respect to AA = Ai — Ao yields a series expansion for 
the difference in chemical potential A/r ea: = ^ X {X\) — fi ex (Xo)\ 


- P^ x = -0u.( 0) [AA (AA + 2A 0 - 2)] + ^^C , „,a o [u«i(0)] . 

71=1 


( 6 ) 


The cumulants C ni A o [Wez(0)] measure the fluctuations of the electrostatic energy u e i( 0) in the state Ao: 


Cl,A 0 [u e l(0)] = (w ei (0)) Ao , 

(7a) 

C 2: X o [Uel{0)} = (Au^)^ , 

(7b) 

C'3,A o [Wez(0)] = (A u ll) Xq , 

(7c) 

C' 4 ,A o [w e z(0)] = (Au^) Ao - 3(AUez)^ , 

(7d) 


where Au e z = u e i( 0) — (u e i{0))\ o . For a Gaussian distribution, all cumulants for n > 3 are zero. If the distribution is not 
Gaussian, the cumulant expansioni4s|4nfinite or higher-order cumulants diverge, i.e., there does not exist an m > 3 such that 
Cm ^ 0 and C n = 0 for all n > m.cjlij Therefore, the use of only mean and variance in eq j(| would be exact if the underlying 
distribution were Gaussian. The self-interaction u B can be seen to correct the mean and the variance of the u e i distributions. 
When identical orders of AA are added in eq^|we obtain corrected cumulants C'. n y. 

C' XM =Ci,a o - 2 (A 0 - 1 )u.( 0 ) 

C' 2Ao =C 2 ,x o -2 Us (0)/P 
C'm,x 0 = Cm,\ 0 for m > 3 . 
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(8a) 

(8b) 

(8c) 





The correction of the variance is independent of the charge state Ao- The correction of the mean on the other hand is 
proportional to the charge factor 1 — Ao and vanishes in the uncharged case Ao = 1. 

We can view eq [] as a Taylor expansion of the chemical potential around Ao, where the cumulants contain the information 
about the derivatives with respect to the coupling parameter A. We can combine this information about the derivatives in 
a x 2 fit °f /r ex (A) to a given functional form. We will here use polynomials pi of varying order l. The constant term, ao, is 
undetermined. The coefficients {a*,} will be chosen to minimize a x 2 functional, 


x 2 (K» = ££[ 4 m) ]" 2 U m) (An,{a fc })-d (m) (A„) 


(9) 


where N is the number of A„ values analyzed. M n is the number of derivatives calculated at A n . p[ Tn ' > is the m-th derivative of 
the polynomial pi( A, {a/c}) = p ex { A) with respect to A. d ^ is the observed derivative, which is related to the cumulant C' m X 
through eqs ^ and |s| 

rf (m) (A) = —/3 m_1 C m ,x ■ (10) 

The estimated statistical error a£ m ' ) of d^ m \Xn) is assumed to be Gaussian. 

We can calculate u e i distributions at different charge states A from simulations of systems with one modified water molecule 
in a solution of unmodified water. For A = 0, the statistical efficiency is greatly enhanced since we can average over all water 
molecules. Using eqs li and the cumulant data of the different charge states can be combined to obtain the chemical 
potential as a function of the coupling parameter A. 


B. Bennett’s overlapping-histogram and acceptance-ratio methods 


We also apply Bennett’s method of overlapping histograms lldl__l to calculate the free energy of uncharging a water molecule. 
We study two states 0 and 1 with a difference in their configurational energies Au = w(Ai) — «(Ao). The normalized probability 
densities Pi(x) ( i = 0,1) of the energy difference Au are defined as 


Pi(x) = (S(Au - x )) A . , 

where 5(x) is Dirac’s delta distribution. 

Canonical-ensemble averages of a function A of the configuration variables in systems 0 and 1 are related through 


(^4)ai — 


(Aexp(- f3Ai)) Xo 


{exp(-(3Au)) Xo 

Using eq |I| we find: 

(A)\ 1 = exp(/3Ap ex ) (Aexp(-/3Au)) ko 
Correspondingly, the probability densities po and pi are related through 

Pi(x) = p 0 (x) exp(0Ap ex - j3x) 


( 11 ) 


( 12 ) 


(13) 


(14) 


with Ap ex = p ex (Xi) — p ex (Xo)- hr a region of overlap between pi and po one expects a linear dependence of \a[p\(x)/po(x)\ 
on x with slope — (3 and intercept /3Afi ex . 

Applied to the case of uncharging a water molecule, we have 


A u = — AX u e i( 0) + A\ (AA + 2Ao — 2) u s (0) . 


(15) 


The probability density po can be calculated from histograms of the electrostatic energy Au = — w e z(0) — m 3 (0) of water 
molecules in bulk water (state 0). To calculate pi, a system has to be studied comprising one uncharged and IV — 1 charged 
water molecules, where u e i is the fictitious electrostatic energy of the uncharged particle obtained by turning on its charges, pi 
is then the probability distribution of Au = —u e i( 0) — m s (0). 1 - 11-1 

From eq M we find the basic relation of Bennett’s acceptance-ratio method to calculate A^i e!E )2H3 


/3Ap ex = In 


{f(-PAu + c)) Xi 
(/(/3Au - c)) Ao 


+ C 


(16) 


where the Fermi function f(x) = 1/[1 + exp(a:)] has been chosen, to minimize the error in the Ap ex estimate, c is an arbitrary 
constant. Setting it to c = 8Ap ex minimizes the expected error .□ In a graphical procedure we search for the intersection of the 
Ao and Ai averages in eq |I^ as functions of c. In addition, we expect a plateau for Ap ex calculated from eq [1^ for c values close 
to the optimum c = /3Ap ex . A more detailed discussion of the two methods due to Bennett can be found in refs [ljand 
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C. Effective electrostatic interactions 


To avoid surface effects, fluid systems are commonly studied using periodic boundary conditions. This results in well-known 
difficulties regarding the treatment of long-range interactions which are large even at distances comparable to the dimen aiana i oj 
the simulation box. We treat the electrostatic interactions in a truly periodic format by using Ewald lattice summation.EJ't3ilJ 
In the Ewald formulation, the electrostatic energy U of a system comprising N water molecules carrying partial charges q a can 
expressed as 


U = 


E EE 9a9 ^ £W ( ri “ J >) + i 

l<i<j<N a=l (3=1 


N 3 3 

EEE qaq/3tpEw{ri a ip) ■ 

i= 1 c=l 0=1 


(17) 


A lattice vector n is added to the distance vector r i a j p such that r i a j /3 = r J/3 — n a + n is in [—L/2, L/ 2] 3 . 

The self-interaction is defined as tpEw(v) = Pew{ r) — l/|r|. The effective Coulomb interaction ifiEW can be expressed in 
terms of rapidly converging lattice sums, 


, , erfc(n|r + n|) 4-7T i k~ \ n 


k#0 


(18) 


V is the volume of the box, erfc is the complementary error function, and k = |k|. The two lattice sums extend over lattice 
vectors n and k of real and Fourier space, respectively. 

To calculate the free energy of charging, we need expressions for the electrostatic energy of the system when a single water 
molecule carries modified charges (1 — X)qo and (1 —A )qH- We will identify the difference of the interaction energies u(Ai) — u(Ao) 
with the difference of the total energies (eq [g]) for the two A values. Correspondingly, we can write for the interaction energy 
m(A) of the modified water molecule i = 1: 


N 3 3 3 3 

U W = E E E^ 1 “ x )9^VEw(ri ajf) )+2 <>-«’££ q a q0ipEw{ri a i p ) ■ (19) 

j =2 ol= 1 (3=1 ot=l (3=1 

We identify the first and second sum with u e i( A) and u s { A), respectively. Eqjl^ contains self-interactions owing to the presence 
of lattice images. We can average the-i»elf4pteractions over the solid angle assuming isotropy. The Ewald potential ipEW can 
be expanded in spherical harmonics£3ii3’E-J In an average of if>Ew{ r) over all orientations, the orthogonality of the spherical 
harmonics yields only a constant plus an r 2 term, 


where the constant is £ew = 


Jn 

—2.837297/L0 From eqs |l| and [lij we obtain for 
«s(A) = -( 1 -A ) 2 ^3 rn 


the self-int« 


( 20 ) 


( 21 ) 


where m is the dipole moment of a water molecule. 

In common implementations of Ewald summation, the Eourier-space sum is rewritten such that it scales with the product 
of the number of particles times the number of k vectors.U These implementations have the disadvantage that electrostatic 
potentials at the charge sites are not directly available^— Ta ciscumvent this problem we use an expansion of the Ewald potential 
in terms of harmonic functions with cubic symmetryH3Ejc3 to calculate the electrostatic potentials (j>o , 4’h 1 , and 4>h 2 - We 
include in the expansion kubic harmonics up to tentE-order. The coefficients are those previously used in the calculation |Of 
chemical potentials of restricted-primitive-model ionsJlj They are listed in Table [E] in the convention of Adams and DubeyO 
For the calculation of the interaction energies in the Monte Carlo simulations, we do not use the kubic-harmonic approximation 
but conventional lattice sums. nn 

We also use a generalized reaction-field (GRF) method for the electrostaticsEjO The GRF interaction depends only on the 
distance and has a cutoff length r c , 

Pgrf{t) = - p(r/r c ) 0(r c - r ) . (22) 

r 

0 is the Heaviside unit-step function; p(x) is a polynomial: 

p(x) = (l-x) 4 (l + 8x/5 + 2x 2 /5) . (23) 

By analogy with Ewald summation (eq §) , we define the A-dependent energy as 
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(24) 


N 3 3 3 3 

=eeE(* - x)qc,q^GRF{ri ajp )+-(i - a) 2 EE gag/3 4'GHF(ri c< l /3 ) . 

3=2 a=l /3=1 Q = 1 (3=1 

The self-interaction is defined as 'kGKF(r) = ipgrf(v) — 1/r. The GRF self-term for an uncharged molecule contains a dominant 
dipolar contribution —4m 2 /r 3 and additional terms of higher order in r^ 1 . 

The direct electrostatic energy u e i is a sum of pair-additive interactions with the other molecules. It closely resembles the 
expression from standard electrostatics, with the interaction 1/r being replaced by ip(r). Based on this resemblance, u e i could 
form a possible choice in the calculation of energies of charging Am. The self-term u s is a consequence of the difference of 
the effective interaction <p(r) and 1/r. In the following, we will consistently include the self-term u s in the calculation of 
electrostatic energies. Taking only u e i would be correct for an isolated system, but not when periodic boundary conditions 
are taken seriously. The direct electrostatic energy u e i stems from the interactions with the other particles in the box (and, 
in the case of Ewald summation, their periodic images). The self-term « S (A) is the difference of the self-interactions with i p as 
effective potential and the vacuum self-interactions with a 1/r potential (which cancel the infinities of the former). In the case 
of Ewald summation, u s has physical significance as it accounts for the interactions with the images of the water molecule with 
charges (1 — A )q a - In the context of GRF electrostatics, u s is introduced by analogy. 

In a previous simulation study, Rick and BerneO also used Ewald summation to calculate the free energy of water but with 
a potential corresponding to vacuum boundary conditions. Their effective potential for u e i calculations differs by —2ivr 2 /3V 
(and a constant) from the Ewald potential used here. Correspondingly, the self-interaction of eq ^ojthen reduces to a constant 
and the self-term of Rick and Berne for neutral molecules is zero. 


III. COMPUTER SIMULATIONS 

We study bulk water (A =p-0) and a system comprising one uncharged and N — 1 unmodified water molecules (A p-fta We 
use the SPC model of water E j Configuration-space averages are performed using the Metropolis Monte Carlo methodEj’Ej The 
systems are studied in the canonical ( NVT) ensemble. The number density is p = 33.33 nm -3 . The temperature is kept at 
T = 298 K. The Monte Carlo move widths is chosen to give an approximate acceptance ratio of 0.5. A cubic box under periodic 
boundary conditions is used for the simulations. The system size is varied between N = 256, 64, and 32 water molecules to 
study finite-size effects. i— ii—ii— i 

The long-range electrostatic interactions are treated using Ewald summation.liil’t3iiJ The real-space screening factor is set to 
rj = 5.6 /L (L is the length of the box). The real-space interactions are truncated at L/ 2.i-The Fourier-space sum is spherically 
truncated using k vectors with |k | 2 < 38(27r/L) 2 , resulting in 2 x 510 k vectors considered.^ The background dielectric constant 
is corrected from infinity to crf = 65 by adding a term 27rM 2 /(2e.RF + 1)V to the potential energy, where M is the total dipole 
moment of the simulation box. 

We also perform simulations using the GRF interaction. A cutoff of r c = 0.9 nm (N = 256) and r c = L/2 (N = 64) is used. 
Again, a correction term for a finite dielectric background is added to the total energy. In all simulations, periodic boundary 
conditions are applied with respect to atomic sites. As starting structures we use randomly oriented and positioned water 
molecules or final structures or previous runs. Before averaging, the systems are extensively equilibrated. 

To calculate the electrostatic energy of charging an uncharged water molecule, systems are studied comprising N — 1 SPC 
water molecules and one uncharged Lennard-Jones sphere with water parameters. The number density is p = 33.33 nm -3 . 
At regular intervals, the electrostatic energy is calculated for 50 random orientations of a fictitious charged molecule with the 
oxygen position identical to the uncharged Lennard-Jones particle. 


IV. RESULTS AND DISCUSSION 

We will first present results for the statistical distribution of u e i. Figure |l| shows histograms of u e i for the charged and 
uncharged state calculated using N = 256 particles and Ewald summation. The histograms are shown on a logarithmic scale so 
that Gaussians appear as parabolas. Also shown are Gaussian distributions with identical mean and variance. The distribution 
in the charged case (A = 0) is approximately Gaussian, but more centered. It decays faster than Gaussian and has a negative 
kurtosis. The distribution in the uncharged case (A = 1) on the other hand deviates more strongly from a Gaussian form. It is 
skewed, has weakly decaying tails and, correspondingly, a positive kurtosis. 

Tabic [n| lists the cumulants C n ,\ of u e i for simulations with A = 0 and 1. Results are shown for Ewald-summation and 
GRF electrostatics and for system sizes of N — 32, 64, and 256 molecules. The results for small system sizes are included 
to study the dependence of the cumulants on electrostatic interaction and system-size, as we are ultimately interested in the 
thermodynamic limit. The consideration of the self-term greatly reduces the variation for the averages Ci,x=o- The uncorrected 
Ewald-summation results for N = 32, 64, and 256 range between —96.69 and —98.00 kj mel -1 . With the self-term added, the 
range is —98.05 to —98.17 kj mol -1 , compared to a typical statistical error of 0.1 kj mol -1 .E3 The GRF self-term is significantly 
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larger than that of Ewald summation. The corrected value C[ A=0 of the GRF simulation with N = 256 agrees with the Ewald- 
summation result within the statistical errors. The corrected TV = 64 result for C[ A=0 is somewhat too negative, but in much 
better agreement with the Ewald-summation data (1 versus 10 kj mol” 1 difference). 

On the other hand, the variance data C^x from Ewald-summation simulations do not show a clear improvement by addition 
of the self-term. However, this is not conclusive since the corrections are of the size of the statistical error or smaller. The 
improvement by addition of the self-term is more evident in the GRF variances. The higher cumulants C 3 ,x=o and C 4 ,x=o show 
less system-size dependence. Results for different system sizes and electrostatics agree within statistical errors. A possible 
exception is that the GRF values for C 3 ,x=o are too small. Overall, the first four cumulants of the charged state (A = 0) can 
be calculated accurately for systems with as few as IV = 64 (and possibly N = 32) water molecules using Ewald summation if 
the finite-size corrections are applied. 

We use the cumulants C[ A , C' 2 \, and C<± } \ for the charged (A = 0) and uncharged (A = 1) case of the N = 256 

simulations to calculate the chemical-potential difference Ap ex of uncharged and charged water. The cumulants C^x for A = 0 
and 1 (see Table |n|) differ significantly. Therefore, we have to use fitting polynomials of order five or higher in the \ 2 fit eq 
Results for Ap ex are listed in Table flli. The polynomial pg of order 8 is interpolating, i.e., pg fits the data exactly. The results 
for polynomials ps to pg vary between Ap ex = 35.13 and 35.80 kj mol” 1 , or by about 2 %. The Ewald and GRF results for like 
orders differ by less than 0.09 kj mol” 1 (0.25 %). The interpolating polynomials pg yield 35.60 (Ewald) and 35.63 kj mol” 1 
(GRF) for Ap ex . From block averages, we estimate the statistical errors of the pg data for Ap ex to be 0.15 kj mol” 1 . n 

We will now compare the results from polynomial fits with those obtained from Bennett’s method of overlapping histograms.u 
Figure shows the ratio ln(pi/po) of the probabilities as a function of A u = u(Ai) — w(Ao), as described above. The distributions 
pi and po are calculated from histograms of u e i using a bin widths of 0.1 kj mol” 1 . Self-interactions are added and the 
sign is inverted according to eq [l5| In the overlap region of the distributions, we indeed find the expected linear behavior 
ln(pi/po) = (3(Afi ex — Ait). We fit constants Ap ex = 35.63 and 35.60 kj mol” 1 to the Ewald-summation and GRF data for 
20 < Ait < 50 kj mol” 1 , respectively. The Ap ex values are in excellent agreement with A^t eiE = 35.60 (Ewald) and 35.63 kj mol” 1 
(GRF) calculated from the polynomials pg interpolating the derivative data. 

Also included in Figure Q are the results of Gaussian approximations to the probability densities. However, as expected 
from the significant deviations from Gaussian behavior found in Figure [j], we do not observe a linear regime. This provides 
evidence for the failure of a simple Gaussian picture of the fluctuation statistics for the hydration of water. Assumption of 
Gaussian behavior (or, correspondingly, quadratic dependence on the coupling parameter) results in inaccurate estimates of 
the chemical-potential difference Ap e: jk 

Bennett’s acceptance-ratio methodn yields 35.59 (Ewald) and 35.60 kj mol” 1 (GRF) for Ap ex . These values are calculated 
by searching for the intersection of the Fermi-function averages for A = 0 and 1 in eq [nj. The corresponding value of c is 
identified with the difference in the chemical potential, /3Ap ex = c. The averages are calculated from the tabulated histogram 
data for An. The acceptance-ratio values of Ap ex are again in excellent agreement with those calculated using the interpolating 
polynomial pg. p, 

Rick and Bernetl calculated Ap ex by evaluating Ci : \ at different charge states —0.22 < A < 1 and combining the data by 
numerical integration. This is essentially a variant of the method we are using in this work. However, inclusion of accurately 
sampled cumulants of higher order—as is done here—should decrease the number of A points required. From eleven molecular 
dynamics simulations of 40 ps duration with N = 512 SPC water molecules, Rick and Berne find Ap ex = 35.15 ± 2.1 kj mol” 1 , 
in agreement with our data. 

Figure shows the chemical potential Ap ex as a function of the charge coupling parameter A. Results are shown for the 
interpolating polynomials p»iOf the Ewald-summation and GRF data, which are practically indistinguishable. Also shown are 
the data of Rick and Berne, u which are found to agree with the pg curves. We find agreement also in the region A < 0, where 
the polynomials pg are extrapolating the data at A = 0 and 1. 

The parabolic Ap ea: curves obtained by assuming Gaussian-fluctuation statistics are included in Figure []. The quadratic 
expansions around both the charged and uncharged state are accurate only for small perturbations | AA| < 0.3 — 0.4. For larger 
perturbations, the parabolas show large deviations from the fitted polynomials and the data of Rick and Berne. This provides 
further evidence for the failure in the large-perturbation regime of simple models, which represent the chemical potential as 
quadratic functions of the coupling parameter based on the assumption of Gaussian-fluctuation statistics. 

Also shown in Figure [| are the results of direct calculations using eq which can be rewritten for the case of uncharging 
water as 


p ex (Ai) - p ex (A 0 ) = -fcsTln (exp[/3(Ai - A 0 )u e i(0)]) Ao + AA(AA + 2A 0 - 2)u s (0) . (25) 

Correspondingly, the direct calculation determines the logarithm of the characteristic function of the distribution of electrostatic 
energies u e i(0). We can see that the exponential average for both Ao = 0 and 1 works well for small perturbations |Ai — Ao| < 0.5. 
For larger perturbations, the direct method fails. Expansions around the charged and uncharged state predict values for the 
difference in chemical potential which are too large by about 2.8 and too small by about 6.8 kj mol” 1 , respectively. 

However, by fitting polynomials to the derivatives at the two ends A = 0 and 1 of the interval, we essentially attempt to 
connect the two curves smoothly near A = 0.5. This explains the success of the method based on fluctuation statistics: We can 
accurately combine information of two widely separated states by smoothly connecting the curves from direct averages near 
A = 0.5, where they are still accurate. Using higher-order cumulants which can be sampled accurately in this case since they 
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are significantly different from zero allows us to construct a non-trivial form of the chemical potential as a function of charge. 
This diarging function deviates strongly from a simple quadratic form, in agreement with previous observations by Rick and 
Bernal and with large values of the kurtosis C 4 /C 2 of the u e i distributions. 

The idea of smoothly connecting the curves from direct calculations can be turned into an algorithm. We shift the A = 1 
curve 71 (A) = — kgT ln(exp[— f3u(\) + /3m(1)])i vertically by Ap ex and determine the intersection with the A = 0 curve 70 (A) = 
— fcBTln(exp[— pu(X) + /3u(0)])o- We then calculate the first derivatives of the two curves at the intersection from polynomial 
interpolation. The absolute value of the difference A = 71 (A) — 70 (A) between the two derivatives is minimized with respect to 
the vertical shift, 

min 1 7 ! (A)- 70 (A) | , (26) 

where A is the solution of 70 (A) = 71 (A) + /Sp BX . |A| is plotted in Figure ^ as a function of the chemical-potential difference 
Ap ex between the charged and uncharged state. The difference in the derivatives shows a distinct minimum for Ap ex between 
35.7 and 35.8 kj mol -1 . This result is in excellent agreement with the previously calculated numbers, deviating by less than 
0.1 %. This simple, geometrical analysis can therefore complement the more elaborate methods and provide accurate results 
for the chemical potential. 

We can also use this analysis to choose an optimal A value in the interval [Ao,Ai] to improve the accuracy of the results. 
From the |A| dependence on the A-value of intersection of the two curves 70 (A) and 71 (A) + A pl x , we find that the smoothest 
connection occurs for A ~ 0.45 — 0.5. For the purpose of calculating free energies, it would be interesting to perform a single 
simulation at A = 0.5 (i.e., for a water molecule carrying 0.5 times the full charges). Judging from the previous results, a 
simulation at A = 0.5 should produce sufficient information to calculate the chemical potential on the whole interval 0 < A < 1. 
However, in this work we restrict our analysis to the physically more relevant systems of bulk water (A = 0) and an uncharged 
solute in water (A = 1). 

Another observation of this study is the importance of the inclusion of self-terms as approximate corrections for finite-size 
effects. Included in Table |ni| are the Ap ex values calculated from polynomial hts to the data of the N = 64 simulations using 
Ewald summation. The values of the N = 64 and N = 256 systems agree closely. The difference of the ps data is 0.04 kj mol -1 
compared to estimated statistical errors of 0.15 kj mol -1 . On the other hand, the ps values for the uncorrected cumulant 
data (i.e., without self-terms) are 35.22 (N = 64) and 35.52 kj mol -1 (N = 256), differing by 0.3 kj mol -1 . Addition of the 
self-term results in a reduced system-size dependence. 

The improvement is even more evident when we compare results from simulations using different electrostatic interactions. 
The difference in chemical potential Ap, ex of charged and uncharged state calculated from the polynomial fit ps and Bennett’s 
overlapping-histogram and acceptance-ratio methods agree within 0.04 kj mol -1 for GRF and Ewald-summation data if self¬ 
terms are added. Without the self-terms they would differ by about 1.7 kj mol -1 . In addition, also the chemical potential as a 
function of the charge shown in Figure is identical in the range of A considered. This shows that even with less sophisticated 
methods for the electrostatics such as GRF compared to Ewald summation it is possible to calculate charge-related free energies 
accurately if only self-terms are considered consistently. 


V. CONCLUSIONS 

Using fluctuation statistics of two equilibrium simulations (bulk water and water hydrating an uncharged Lennard-Jones 
particle) we are able to calculate accurately the chemical potential of SPC water as a function of its charge. Bennett’s 
overlapping-histogram and acceptance-ratio methodsEI and a geometrical method based on smoothly connecting the curves of 
direct, exponential averages at the two states give identical results-for the chemical-potential difference between charged and 
uncharged water. Our results agree with those of Rick and Bernctl calculated from thermodynamic integration using eleven 
states with respect to both the chemical-potential difference and the dependence on the charge. 

We find signfficant deviations from a quadratic dependence of the chemical potential on the charge coupling parameter. This 
has important implications. It shows that even for the fairly simple system of water in water second-order perturbation or, 
equivalently, assumption of Gaussian-fluctuation statistics, allow only crude descrip4jeuis-|Of the actual thermodynamics. This 
also affects the potential usefulness of linear continuum models of electrostatics, □oEJlj which by design do not go beyond 
a quadratic behavior. However, when calculating free-energy differences, continuum models usually compare approximate 
representations of physical states rather than performing an expansion around a single state. This can explain their success 
of giving at least approximately correct free-energy values even in the presence of strongly non-linear dependencies in a 
corresponding system of atomic resolution .□ 

In addition to thcfrailure of the simple Gaussian model, expansions around only a single state using higher-order cumulants 
are expected to failEa Increasing the order of the perturbation requires accurate information about the poorly sampled tails of 
distributions. This merely reflects the difficulty to extrapolate to states that strongly differ in their structure and fluctuation 
statistics. Interpolation on the other hand is in general a much simpler task and allows more accurate predictions. In this work 
we effectively combine the information of two states to derive a polynomial expression for the chemical potential. 
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An important point concerns the treatment of electrostatic interactions in computer simulations of systems under periodic 
boundary conditions. With respect to system-size dependence and electrostatic model, we obtain consistent results using Ewald 
summation and a generalized reaction-field model. But consistency is only achieved if self-interactions are included. These 
self-interactions arise naturally when effective potentiaknare used for the Coulomb interactions. We made similar observations 
in previous studies of the chemical potential of ionsotj Adding self-terms to the energies acts as an effective correction for 
effects of a finite system size. Neglecting them can result in significant deviations for small systems of a few hundred particles. 
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FIG. 1. Probability density of the electrostatic energy u e i of a water molecule with full charge (A = 0; centered near 
uei = —100 kj mol -1 ) and with zero charge, where fictitious charges are turned on (A = 1; centered near u e i = 0). The 
probabilities are shown on a logarithmic scale such that Gaussians appear as parabolas. Also shown are Gaussian distributions 
(dot-dashed lines) with mean and variance equal to the calculated distributions. The distributions are calculated from Monte 
Carlo simulations using N = 256 particles with Ewald summation. The curves do not contain corrections for self-interactions. 


FIG. 2. Bennett’s method of overlapping histograms.u The ratio ln(pi/po) of the probabilities is shown with crosses as a 
function of A u = u (Ai = 1) — u(Ao = 0). pi and po are approximated by histogram values calculated from Ewald-summation 
and GRF simulations using N = 256 particles. The GRF data are shifted vertically by —10. Also shown are lines with slope 
—(5 fitting the data. The dot-dashed lines represent the results of a Gaussian approximation (see Figure [lj) with mean and 
variance taken from the GRF and Ewald-summation data. 


FIG. 3. The chemical potential Ap ex of water as a function of its charge. A = 0 and 1 correspond to the fully charged and 
uncharged state, respectively. Solid lines are polynomials of order 8 fitted to the cumulants C[ x , C ' 2 63 ,a, and C 4 ,a for A = 0 

and 1. The curves for Ewald-summation and GRF data are practically indistinguishable. Shown as symbols with error bars 
are the data of Rick and Berne taken from Figure 2 of ref Also shown are results from direct calculations using eq ||^. The 
expansion around charged and uncharged state A = 0 and A = 1 are shown with dashed and dot-dashed lines, respectively, in 
both cases using the data of Ewald-summation simulations. The dotted lines represent quadratic expansions around A = 0 and 
1 using mean and variance of the Ewald-summation data and assuming Gaussian statistics. 


FIG. 4. Smooth connection of the chemical-potential curves. Two curves for the chemical potential are calculated from direct, 
exponential averages eq for the charged and uncharged state. One of the curves is shifted vertically and the intersection 
of the two curves is calculated. Shown is the absolute value of the difference A of the first derivatives at the intersection 
as a function of the corresponding difference in chemical potential Ap ex . Where |A| reaches a minimum, the two curves can 
be connected most smoothly. This determines an estimated value for the chemical-potential difference between charged and 
uncharged water Afj, ex = 35.7 — 35.8 kj mol -1 . 


TABLE I. Expansion-coefficient*. S and A; of the Ewald potential pew{ r) in terms of kubic-harmonic functional!} in the 
convention of Adams and Dubey.Efl The constant S has been adjusted from its exact value such that the numerical scheme 
gives a vanishing average potential of a point charge in a cube. 


5 

-2.8373002368 

A 2 

2tt/3 

A 4 

7.718196 

Aq 

20.657378665 

As 

86.85346475 

A 10 

179.631024892 


TABLE II. Cumulants of the u e i distributions calculated from Monte Carlo simulations using N particles. The cumulants 
were averaged over P passes, where one pass consists of N attempted moves (i.e., one attempted Monte Carlo move per particle). 
Coulomb denotes the treatment of electrostatic interactions. The fully charged and uncharged water molecules correspond to 


A = 0 and A = 1. C k denotes the fc-th cumulant of the distribution [measured in (kj mol 1 ) fc ], The cumulants with corrections 
for the self-interaction are listed as C' k . 

N 

Coulomb 

P/1000 

A 

Gi 

c\ 

c 2 

C ' 2 

c 3 

c 4 

256 

EW 

280 

0 

-98.00(10) 

-98.17(10) 

439.5(2.0) 

440.0(2.0) 

120(50) 

-25400(1000) 

64 

EW 

860 

0 

-97.47(10) 

-98.15(10) 

439.7(2.0) 

441.4(2.0) 

102(60) 

-24900(1500) 

32 

EW 

860 

0 

-96.69(15) 

-98.05(15) 

443.4(3.0) 

446.8(3.0) 

100(90) 

-27200(2500) 

256 

GRF 

460 

0 

-94.73(10) 

-98.22(10) 

434.5(1.5) 

443.1(1.5) 

62(50) 

-24400(1000) 

64 

GRF 

900 

0 

-88.36(10) 

-99.05(10) 

418.3(1.5) 

444.7(1.5) 

21(40) 

-23650(1000) 

256 

EW 

280 

1 

-0.006(5) 

-0.006(5) 

110 . 0 ( 2 . 0 ) 

110.4(2.0) 

261(20) 

13000(1000) 

64 

EW 

1000 

1 

-0.007(5) 

-0.007(5) 

109.5(3.0) 

111.2(3.0) 

263(15) 

12500(1000) 

256 

GRF 

600 

1 

0.16(2) 

0.16(2) 

106.4(2.0) 

115.1(2.0) 

277(20) 

12450(1000) 
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TABLE III. Difference in chemical potential Ap ex of an uncharged and charged water molecule in water. Results from x 2 
fits of polynomials p n of order n to the derivative data of the N = 256 simulations using Ewald-summation (EW-256) and 
GRF interactions. Also included are results for the N = 64 simulations using Ewald summation (EW-64). A^i ex is listed in 
kj mol -1 . The statistical errors of the p$ data are estimated to be 0.15 k.I mol. 



EW-256 

EW-64 

GRF 

P5 

35.46 

35.23 

35.45 

P6 

35.13 

35.05 

35.22 

P7 

35.76 

35.71 

35.80 

PS 

35.60 

35.56 

35.63 
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Fig. 2 (Hummer, Pratt, and Garcia) 
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Fig. 4 (Hummer, Pratt, and Garcia) 



